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-xf : 

The WHOT-QCD Collaboration is pushing forward a series of lattice studies of QCD 
at finite temperatures and densities using improved Wilson quarks. Because Wilson- type 
w ' quarks require more computational resources than the more widely adopted staggered-type 

quarks, various theoretical and computational techniques have to be developed and applied. 
I , In this paper, we introduce the fixed-scale approach armed with the T-integration method, 

|- the Gaussian method based on the cumulant expansion, and the histogram method combined 

with the reweighting technique. Adopting these methods, we have carried out the first study 
i-C ' of finite-density QCD with Wilson-type quarks and the first calculation of the equation of 

state with 2+1 flavors of Wilson- type quarks. We present results of these studies and discuss 
j ■ perspectives towards a clarification of the properties of 2 + 1 flavor QCD at the physical point, 

' at finite temperatures and densities. 

t> : 

■<^j- ■ §1. Introduction 

m : 

Temperature T and density are controllable parameters of the system. At suffi- 
ciently high T, we expect that the confinement is violated and the chiral symmetry is 
. recovered because the effective coupling at the thermal energy scale becomes small 

y— I \ due to the asymptotic freedom of QCD. Then, the systems with quarks and glu- 

j> ■ ons will form a colored plasma state called "quark-gluon plasma" (QGP). Although 

the humankind has never experienced the QGP, QGP is expected to play an im- 
portant role in the creation of matter during the early development of the Universe. 
Furthermore, QGP is considered to be observed by relativistic heavy-ion collision ex- 
periments at RHIC and LHC. 1 ) Similar to the case of high temperatures, we expect 
deconfinement at sufficiently high density too because the average distance between 
quarks becomes small there and the property of the system will be dominated by the 
asymptotic freedom. The density is controlled by the chemical potential [i. At very 
large n and low T, we expect a BCS-like state called "color superconductor" due to 
an attractive interaction among quarks. At lower densities around the nuclear den- 
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Fig. 1. Prospected phase diagram of QCD at finite temperatures and densities. 

sity, we expect a nuclear fluid state, which may appear around the core of neutron 
stars. We thus expect a rich phase structure in QCD as a function of T and \i. See, 
e.g., Ref. 2) for a recent review. A prospected phase diagram is shown in Fig. 1. 

When we vary the quark masses off the physical point, the nature of the quark 
matter may be different depending on them. The usual expectation for the order of 
the finite-temperature QCD transition at fj, = 0, based on effective model studies and 
lattice simulations, is summarized in Fig. 2. Details of the phase diagram as well as 
the nature of the quark matter at finite T and fi are, however, not well clarified yet. 
Because the issue is essentially non-perturbative, numerical studies on the lattice is 
the only systematic way to investigate the phase structure directly from the first 
principles of QCD. 3 )' 4 ) 

Most lattice studies of hot/dense QCD have been done with computationally 
less demanding staggered-type lattice quarks. 3 )' 4 ) In particular, in the study of the 
equation of state (EOS), extrapolations to the physical point and to the continuum 
limit have been achieved only with staggered-type quarks. However, the theoretical 
basis such as locality and universality are not well established with them. 3 ) There- 
fore, to evaluate the effects of lattice artifacts and thus to obtain reliable predictions 
to be compared with experiment, it is important to perform simulations using the- 
oretically sound lattice quarks, such as the Wilson-type quarks. Here we note that, 
until recently, the 0(4) scaling property expected around the chiral transition of two- 
flavor QCD 5 ) has been observed only with Wilson-type quarks. 6 )' 7 ) Quite recently, 
an O(N) scaling behavior*) was observed with an improved staggered quarks by let- 
ting the light quark mass much lighter than the physical u and d quark masses. 8 ) 
Therefore, some of the chiral properties around the transition temperature may be 



Presumably the 0(2) scaling, as expected from the symmetry of staggered-type quarks. 
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Fig. 2. Prospected order of the finite-temperature QCD transition at /i = as a function of the 

def 

light quark mass m U d =' m u = mj and the strange quark mass m 3 . The top-right corner 
corresponds to the quenched limit of QCD. Lattice simulations with improved staggered quarks 
suggest that the physical point is located in the crossover region. The phase diagram is not fully 
established yet. See Ref. 3) for discussions and caveats. 



easier to extract with Wilson-type quarks. 

A reason that Wilson-type quarks have not been widely adopted in the study of 
hot/density QCD is that the computational cost is larger than that for staggered- 
type quarks, in particular at small quark masses. Previous studies with Wilson-type 
quarks were limited to large quark masses and the case of two-flavor QCD at van- 
ishing chemical potentials. 7 )' 9 ) The WHOT-QCD Collaboration is pushing forward 
studies of lattice QCD at finite T and \x adopting improved Wilson quarks 10 ) coupled 
to RG-improved Iwasaki glues. 11 ) We want to extend the studies to more realistic 
2 + 1 flavor QCD at finite chemical potentials with physical light quarks. Towards 
this goal, we made a series of simulations by implementing and developing efficient 
methods for Wilson-type quarks. We developed the T-integration method to make 
the fixed-scale approach applicable, 12 ) tested the Gaussian method to tame the sign 
problem, 13 ) and extended the histogram method by combining with the reweight- 
ing technique to investigate the phase structure. 14 )~ 16 ) With these techniques, we 
have studied thermodynamic properties of the quark matter through the equation 
of state, 17 )' 18 ) heavy-quark free energies and screening masses, 19 )~ 21 ) spectral func- 



In Sect. 2, we review techniques adopted in our investigation of finite-temperature 



tions, 22 )' 23 ) etc. 
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QCD on the lattice. We also introduce the T-integration method to calculate the 
EOS in the fixed-scale approach. We then report our calculation of EOS with 2+1 
flavors of improved Wilson quarks in Sect. 3. Sections 4, 5 and 6 are devoted to our 
study of finite-density QCD. We first discuss in Sect. 4 major methods to perform 
simulations of QCD at finite densities. We introduce our approach using the cu- 
mulant expansion in a hybrid method of Taylor and reweighting methods. We then 
present in Sect. 5 our results of the pressure and the quark number susceptibility at 
finite densities in two-flavor QCD, adopting these methods. In Sect. 6, we introduce 
another method, a histogram method, to investigate the first order transition and 
its boundary. We apply the method to calculate the location of the critical point 
where the first-order deconfining transition of QCD in the heavy quark limit turns 
into a crossover as the quark masses are decreased. We then study the critical point 
at non-zero chemical potentials. We also present our on-going project to study finite 
density QCD with light dynamical quarks by combining the histogram method with 
phase-quenched simulations. Finally, in Sect. 7, our results of the heavy-quark free 
energy are summarized for zero and finite densities. A short summary is given in 



On a lattice with a size iV s 3 x N t = N s \t e , the temperature of the system is given 
by T = l/(N t a), where a is the lattice spacing. To vary T, we may either vary a at 
fixed Nt, or vary Nt at a fixed a. Let us call the former as the fixed-N t approach, 
and the latter as the fixed-scale approach. 

In the fixed-iVt approach, we can vary T continuously through a continuous 
variation of a. This is a reason that the fixed- Nt approach has been widely adopted 
in many simulations. 

The value of a is controlled by the coupling parameters, which we denote as b. 
For QCD with Wilson-type quarks, b = (/?, k u , Kd, k s , ■ ■ ■ ). We first define the lines 
of constant physics (LCP's) in the coupling parameter space by the fixed dimension- 
less ratios of physical observables such as m^/nip = m^a/rripa. Here, to remove 
additional dependence on T, these observables have to be measured at T = 0. A 
LCP represents a physical system at different values of a. In two-flavor QCD with 
improved Wilson quarks, LCP's defined by m^/nip are determined in Refs. 9) and 
19). In 2+1 flavor QCD, we have to fix one more dimension-less ratio such as m^/m,, 
or rrir^Jm^. Our world corresponds to the LCP with m^/rrip ps 135/770, etc. 

The beta function adb/da is defined as the variation of b along a LCP. In the 
fixed-iVt approach, we vary T of a given physical system by varying b along a LCP 
on a lattice with a fixed value of Nt- 

The energy density e and the pressure p of the system are given by derivatives 
of the partition function Z in terms of T and the physical volume V = (N s a) 3 : 



where (■ ■ ■ ) su b is the thermal average with zero temperature contribution subtracted 



Sect. 8. 



§2. Thermodynamics of QCD on the lattice 




(2-1) 
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for renormalization. To vary T and V independently, we need to introduce anisotropic 
lattices. When a s and a t are the lattice spacings in spatial and temporal directions, 
V and T are given by V = (N s a s ) 3 and T = l/(N t a t ). Then, in principle, (2-1) 
can be evaluated by independent variations of a s and at- However, this requires a 
systematic study of a set of physical observables on anisotropic lattices with varying 
both a s and at, which is, however, quite demanding. 
Here, we note that the combination 



d d 

is given in terms of a uniform rescaling at — h a s — — , which can be realized without 

oat das 

introducing anisotropic lattices. We thus obtain on isotropic lattices 

T I d\nZ \ _T db / dhxZ \ T db / dS\ 



f -~ 3p =~y \«— L= -v"Ta \^rL = v a T a • \aL • {2 ' 3> 

where S is the lattice action. The coefficient adb/da is just the beta-function of 
the system, whose non-perturbative values can be determined by simulations on 
isotropic lattices. Eq. (2-3) enables us to study e — 3p non-perturbatively without 
introducing anisotropic lattices. The combination e — 3p is nothing but the trace of 
the energy-momentum tensor, and called the trace anomaly, e — 3p vanishes for free 
gasses, but will have non-trivial values with interacting matters. 

2.1. Fixed-Nt approach and integration method 

In order to obtain e and p separately, we need one more independent input. 
The most widely adopted is the integration method, 2 ^ with which we can determine 
the pressure p non-perturbatively through an integration in the coupling parameter 
space: 

T f\ r I 1 3Z\ T f 6 -. /ds\ 

This relation is obtained by differentiating and then integrating the thermodynamic 
relation p = (T/V) In Z in the coupling parameter space. The integration path is 
free to choose as far as the initial point &o is located in the low temperature phase 
such that p(bo) ~ 0. See Appendix A of ref. 9) for a concrete demonstration of the 
path-independence. 

Several points to be kept in mind with the fixed- Nt approach are as follows: (i) 
When we fix N s , the spatial volume V = (N s a) 3 is varied simultaneously as we vary 
T. In the high T region, V may be quite small with a fixed N s . To keep V around a 
fixed value, N s has to be increased as we increase T. (ii) At low T's, the lattice may 
be coarse. To ensure asymptotic scaling around the QCD transition temperature, a 
large N t together with improvement of the lattice action is mandatory, (iii) For the 
zero-temperature subtraction, we have to carry out zero-temperature simulations at 
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all of the finite-temperature simulation points. Together with the systematic zero- 
temperature simulations to determine LCP in a wide range of the coupling parameter 
space, an indispensable fraction of the total computational cost is required to carry 
out the systematic zero-temperature simulations in the fixed- N t approach. 

2.2. Fixed-scale approach and T -integration method 

In the fixed-scale approach, T is varied by varying Nt at a fixed b (and thus at 
a fixed a), i.e., the simulations are done at the same b point for all T's. Therefore, 
all the simulations are automatically on a LCP without fine-tuning. Furthermore, 
the T = subtractions can be done by a common T = lattice. We may even 
borrow high statistic configurations at T = on the International Lattice Data 
Grid (ILDG). 25 ) Therefore, we can largely reduce the cost for the zero-temperature 
simulations with the fixed-scale approach. 12 ) 

On the other hand, the conventional integral method to obtain p by an inte- 
gration in the coupling parameter space is inapplicable, because data are available 
at only one b point in the fixed-scale approach. To overcome the problem, we have 
developed a new method, the T -integration method: 12 ^ Using a thermodynamic re- 
lation at vanishing chemical potential, we obtain 

T A fM _ i^p i_ _ f T dT i^i (2 . 5) 

8T \T 4 ) ~ T 4 T 4 ~ J To T 5 1 ' 

with p(Tq) ~ 0. When we vary T by varying b along a LCP, the integral in (2-5) 
is equivalent to that in (2-4) with the integration path chosen to be on this LCP. 
However, (2-5) allows us to integrate over T without varying b. 

In the fixed-scale approach, various T's are achieved by varying N t . Because 
N t is discrete, we have to interpolate the data with respect to T to carry out the 
integration of (2-5). We need to check the magnitude of systematic errors from the 
interpolation of the trace anomaly in T. Application of the method to finite [i is 
straightforward when we reweight from /j, = 0. 

We find that the fixed-scale approach is complemental to the conventional fixed- 
N t approach in many respects: At very high temperatures, typically at T > l/3a, the 
fixed-scale approach suffers from lattice artifacts as discussed below, while the fixed- 
N t approach can keep Nt finite and can keep the lattice artifact small by adopting 
a sufficiently large N t . On the other hand, at small T, typically at T < T c , the 
fixed-scale approach can keep a small at the price of larger cost due to large Nt, 
while the fixed- N t approach suffers from lattice artifacts due to large a. 

Another attractive point of the fixed-scale approach in a study with improved 
Wilson quarks is that, unlike the case of the fixed-A^ t approach, we can keep the 
lattice spacing small at all temperatures and thus can avoid extrapolating the non- 
perturbative clover coefficient csw to coarse lattices on which the improvement pro- 
gram is not quite justified. 

It should be kept in mind that the fixed-scale approach is not applicable at very 
high temperatures where, besides the artifacts due to a quite small value of N t , the 
lattice spacing a may also be too coarse to resolve thermal fluctuations. 12 )' 21 ) To 
get an idea about the latter effects, we estimate a typical length scale of thermal 
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Fig. 3. Test of the fixed-scale approach armed with the T-integration method in quenched QCD. 12 - 1 
Left trace anomaly on an anisotropic lattice (a2) compared with the isotropic lattice with 
similar spatial lattice spacing and volume (i2). Right: energy density and pressure by the T- 
integral method. The shaded curves represent the results of the conventional fixed- iVt method 
at Nt = 8. 26) 



fluctuations by the thermal wave length A ~ 1/E where E is an average energy of 
massless particles at finite T. We find E ~ 3T£(4)/£(3) ~ 2.7T for the Bose-Einstein 
distribution and E ~ 3T£(4)/£(3) x 7/6 ~ 3.15T for the Fermi-Dirac distribution. 
We then obtain A ~ 1/3T. Thus, data at T > l/3a for which a > A should be taken 
with care. 21 ) 

We have tested the fixed-scale approach and the T-integration method in quenched 
QCD. 12 ) Main results are summarized in Fig. 3. Comparing EOS' obtained on vari- 
ous lattices as well as the result from the fixed- Nt approach on large lattices, we find 
that the fixed-scale approach is reliable and powerful to calculate EOS, in particular 
at low and intermediate temperatures. The systematic errors due to the interpolation 
in T is well under control in these studies. The EOS from the fixed-scale approach 
was shown to be well consistent with that from the fixed- Nt approach with large 
At (At 8), except for the high temperature limit where the fixed-scale approach 
suffers from lattice discretization errors. 

We adopt the fixed-scale approach to calculate the EOS in 2 + 1 flavor QCD in 
Sect. 3. We also compute heavy-quark free energies in 2 + 1 flavor QCD with the 
fixed-scale approach in Sect. 7.2. 

§3. Equation of state in 2 + 1 flavor QCD with improved Wilson quarks 

A systematic study of finite temperature QCD with improved Wilson quarks 
was made by the CP-PACS Collaboration around the beginning of this century 
for the case of two-flavor QCD at vanishing chemical potentials using the fixed- At 
approach. 7 )' 9 ) From a series of systematic simulations, they determined the phase 
structure and LCP's, confirmed the 0(4) scaling, and obtained the EOS along several 
LCP's in the range m n /m p > 0.65 around and above the pseudocritical temperature 



8 



S. Ejiri, K. Kanaya and T. Umeda 



T pc on N t = 4 and 6 lattices. 

We extend the study to 2 + 1 flavor QCD. By adopting the the fixed-scale ap- 
proach, we use zero-temperature configurations generated by the CP-PACS+JLQCD 
Collaborations. 27 )~ 29 ) Our lattice action consists of the RG-improved Iwasaki gauge 
action 11 ^ and the clover-improved Wilson quark action: 10 ) 

S 9 + S q (3-1) 

-p E { E + E c * w *% } . (3-2) 

x \fi>v n,i> ) 

X Y^uiMid- (3-3) 

f=u,d,s x,y 

S *,y ~ K f E {(* " u ^ 5 x+d,y + (1 + 7m) u l-w 5 x-ii,v} 
- 5x,j/ csw(/3) «/ ^ cr^F^ (3-4) 

where ci = —0.331 and Co = 1 — 8ci. We set n u = Kd = n u d and adopt the clover 
coefficient csw nonperturbatively determined by the Schrodinger functional method 
in 29). F x ^ u = (f x ,ui/ — fl,fii/) /8i is the lattice field strength, with f x ^ u the standard 
clover-shaped combination of gauge links. Hadronic properties with this action have 
been studied down to the physical point by the CP-PACS, JLQCD and PACS-CS 
Collaborations. 27 ^ 33 ) 

As the first determination of the EOS with Wilson- type quarks in 2+1 flavor 
QCD, we study at (/?, k U( i, k s ) = (2.05, 0.1356, 0.1351), that corresponds to the small- 
est lattice spacing and the lightest u and d quark masses (m^/nip ~ 0.63) among 
the zero-temperature configurations generated by the CP-PACS+JLQCD Collabo- 
rations. 27 )' 28 ) The s quark mass is set around its physical point (Tn„ ss /m ( j > ~ 0.74). 
These u and d quark masses are much larger than their physical values yet. A study 
at the physical point 33 ) is reserved for the next step. The hadronic radius at this 
simulation point is estimated to be r^/a = 7. 06(3). 34 ) Setting the lattice scale by 
r = 0.5 fm, we estimate 1/a ~ 2.78 GeV (a ~ 0.07fm). The lattice size is 28 3 x 56 
with N s a ~ 2 fm. 

3.1. Beta functions 

Using the same values of the coupling parameters as the zero-temperature simu- 
lation, we have generated finite temperature configurations on 32 3 x N t lattices with 
N t = 4, 6, ■ ■ ■ 16. 17 )' 18 ) This range of N t corresponds to T ~ 170-700 MeV. 

To evaluate the trace anomaly using (2-3), we need the beta functions a{d(3/da) 
and a(dnf/da) (/ = ud and s). These beta functions can be determined nonper- 
turbatively through the coupling parameter dependence of zero-temperature observ- 
ables. We use the data of am p , m^/nip and m r]ss /m ( j ) at 30 simulation points of 
the CP-PACS+JLQCD zero-temperature configurations 27 )' 28 ) to extract the three 
beta functions. The first observable am p sets the scale. A naive method to ob- 
tain the beta functions is to fit the data of these observables as functions of the 
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Fig. 4. Determination of the beta functions in 2+1 flavor QCD. 18 - 1 Left: The global fit for /9 as 
a function of m p a with corresponding m p /m w and m r)sB /m ( f ) . Square symbols show coupling 
parameters in the CP-PACS+JLQCD study. To avoid a too busy plot, only half of the data 
points are shown. Right: Beta functions on our LCP, m^/mp — 0.6337 and m I)ss /m</ ) = 0.7377, 
as functions of j3. The scale setting is made with am p . Beta functions for K u a an d k s are 
magnified by factor 100. 



coupling parameters (/3, K u d, k s ), and invert the matrix of the slopes d(am p )/d/3 
etc. However, because a{dnt /da)'s are numerically much smaller than a{dj3/da), the 
formers get large relative errors through the matrix inversion procedure, i.e., errors 
for large components contaminate and dominate the errors for small components. 
On the other hand, from the previous experience of two-flavor QCD with improved 
Wilson quarks in the fixed-iVt approach, 9 ) we expect that, although a(dnf /da)'s are 
small, the quark contribution is important in the trace anomaly. Therefore, a precise 
determination of a{dnf /da)'s is required. 

To avoid the matrix inversion procedure, we instead fit the coupling parameters 
b = ((3, K u d, k s ) as functions of three observables am p , m^/nip = X and m r)as /m ( j ) = 
Y. Consulting the overall quality of the fits, we adopt the following third order 
polynomial function of the observables in this study: 

b = c + c\ (am p ) + c 2 (am p ) 2 + c 3 X + &4X 2 + c 5 (am p )X + cqY + c 7 Y 2 
+ c 8 (am p )Y + c 9 XY + c 10 (am p f + c n X 3 + c 12 Y 3 + c 13 (am p )X 2 
+ c M (am p ) 2 X + C15 (am p )Y 2 + c 16 (am p ) 2 Y + c 17 XY 2 + c 18 X 2 Y 
+ C19 (am p )XY. (3-5) 

The global fits for each component of b are independent and have dof = 10. We find 
reasonable x 2 /dof of 0(1)*). The result for /3 is shown in the left panel of Fig. 4. 

In this study, we define LCP's by m^/mp and m r]as lm$ at T = 0. Therefore, 
the beta functions are extracted as ad/3/da = (am p ) df3/d(am p ) etc., in terms of 

^ We find xVdof = 1.63, 1.08, and 1.69 for the fit of /?, /t ut j, and k s , respectively. Here, \ 2 ls 
evaluated using a standard deviation of each coupling parameter estimated by the error propagation 
rule using the errors of the observables and the partial derivatives of the resulting fitting function 
(3-5) with respect to the observables, neglecting the covariance among the observables. 
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Fig. 5. Trace anomaly (e — 3p)/T 4 , energy density e/T 4 , and pressure 3p/T 4 in 2+1 flavor QCD. 18) 
The thin and thick vertical bars represent statistic and systematic errors, respectively. The 
curves are drawn by the Akima spline interpolation. 



the coefficients ci, c*2, C5, eg, • • • in (3-5). The resulting beta functions for our LCP 
(m^/rrip = 0.6337, m^/m^ = 0.7377) are shown in the right panel of Fig. 4 as 
functions of (3. Their values at j3 = 2.05 are used to determine the trace anomaly. 
As the variable to set the scale, we may alternatively adopt am n , arriK or arriK* 
instead of am p in (3-5). We use this freedom to estimate a systematic error. Taking 
the results from am p as the central value, we obtain 

= -0.279(24)0, = 0.00123(41)0, = 0.00046 (26) O 

(3-6) 

at our simulation point, where the first brackets are for statistic errors, and the 
second brackets are for the systematic error. 18 ) 

3.2. Equation of state 

We now calculate the trace anomaly by (2-3): 

e - 3 p _ N? f d(3/dS\ dK ud I OS \ d Ks /8S\ \ 

T* ~ Ni\ a da\dfj/ suh + a da \d Kud / snh +a da\d Ks } > ^ 



s I sub 



with 



d 

~ ~d(3 




+ E «/ ( E tr^^J^ (Mf)~i ) , (3-8) 



f=u,d,s \x,fx>v 
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c^a^F x ^(M f )-y) I, (3-9) 




' sub/ 

where Cf = 2 for f = ud and 1 for / = s. We evaluate the spatial traces in (3-8) 
and (3-9) by the random noise method with U(l) random numbers. 13 ) The number 
of noise is 1 for each of the color and spinor indices. 

Our result for the trace anomaly (e — 3p)/T 4 is shown by a red thick curve in 
Fig. 5. We note that the peak height of about 7 at T = 199 MeV (N t = 14) is 
roughly consistent with recent results of highly improved staggered quarks obtained 
with the fixed N t approach at N t = 6-12. 35 )' 36 ) The shape of (e - 3p)/T 4 suggests 
that T pc locates between 174 and 199 MeV. 

Carrying out the T-integration (2-5) of the trace anomaly, we obtain the pressure 
p shown in Fig. 5. The energy density e is calculated from p and e — 3p. 

Besides the larger errors, our EOS looks roughly consistent with recent results 
with highly improved staggered quarks near the physical point. 35 )' 36 ) We note that 
our peak is slightly higher. This is consistent with the fact that our light quark masses 
are heavier than their physical values: the experience with improved staggered quarks 
suggests that the peak becomes slightly higher as the light quark masses are increased 
[see, e.g., Ref. 35)]. 

Our final goal is to extend the study towards the physical point, adopting the on- 
the-physical-point configurations generated by the PACS-CS Collaboration. 33 ) From 
the present study at m u d heavier than the physical point, we encountered a couple 
of issues: Errors of EOS in Fig. 5 are larger than those obtained with the fixed- N t 
approach. 35 ^' 36 ) Besides smaller statistics, this is due to the large statistical error in 
(e - 3p)/T 4 at T < 200 MeV, which is caused by the enhancement factor iV t 4 in (3-7) 
(note that S is proportional to Nt). To obtain accurate EOS at low temperatures, we 
need a large statistics of 0(N^)*\ This is, however, an unavoidable step to suppress 
discretization errors, and the issue is common with the fixed- iVt approach. Another 
source of systematic errors in EOS is the limited resolution in T due to the discrete 
variation of N t in the fixed-scale approach. Note that the lattice spacings in full 
QCD studies are usually coarser than those in quenched studies. Furthermore, in 
the present study, N t is limited to be even due to the simulation program set we 
have adopted. To improve the resolution in T, we need simulations at odd values of 
N t and/or a finer lattice spacing a. An alternative way will be to combine results 
at different a's, since we can choose a's fine with the fixed-scale approach and thus, 
after confirming that the discretization effects are sufficiently small in the observables 
under study, we may combine the results at different a's to more smoothly interpolate 
in T. We leave these trials to a forthcoming study with lighter quarks. 



A power of N t is reduced by the averaging over the lattice sites. 
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§4. Finite density QCD on the lattice 

Next, let us move on to the issues of finite density QCD. We consider the action 
given by (3T)-(3-3) with the quark matrix M[ y replaced by 



M ly = Sx,V ~ K f Yl ~ Ux^x+^y + (1 + 7/i) U^Jx-^y } 

{ e "/°(i - 74) u xA s x+ly + e"^«(i + 74 ) r; ; ,,,} 



(4-1) 

Here, the theory at [if 7^ is known to have a serious problem: In a Monte Carlo 
simulation, we generate configurations of link variables {U x ,n} with the probability 
in proportion to the Boltzmann weight (J}/ det M?) e~ Sg . The expectation value of 
an operator 0[U] is then evaluated as an average of O over the configurations, 

iVconf - {U^} 

where k = (k u , kj, ■ ■ ■ ) and /I = (/x u , [id, • • • )■ At /2 = 0, the quark determinants are 
real due to the 75-hermiticity of the quark matrices, (M^)t = 75 M^ 7 5. However, 
when [if ^ 0, the 75-hermiticity relation changes to 

M'( M/ )] =75M'(-/x / ) 75 , (4-3) 

and thus detM-^ becomes complex unless [if = 0. Because the Boltzmann weight 
has to be real and positive, we cannot perform a Monte Carlo simulation at [if 7^ 
directly. 

Various methods have been proposed to study finite density QCD avoiding the 
complex weight problem. However, at present, all of them are applicable in small 
[if /T regions only. In the following subsections, we introduce these methods, mainly 
focusing on the Taylor expansion and reweighting methods we adopt. 

4.1. Taylor expansion method 

The simplest approach to study finite density QCD avoiding the complex weight 
problem is the Taylor expansion method, in which physical quantities are expended 
in terms of [if/T around [1 = 0. 37 ^ 39 ^ For example, the pressure p = (T/V) \nZ is 
expanded as 

p = £ <**m (£)'(£)'(£)'. (4-4) 



rp4 

i,j,k=0 



1 1 d i+j+k \nZ 
Cl ' jk ~ i\j\k\ VT 3 d([i u /TYd([i d /T)id{[i s /TY 



(4-5) 

ft=Q 
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in three-flavor QCD. The Taylor coefficients Cijt can be evaluated by a conventional 
Monte Carlo simulation at jl = which is free from the complex weight problem. 
We expect that QCD in the high temperature limit is well described by a free gas 
of quarks and gluons, in which p consists of terms proportional to /j/j and /it 4 only. 
Therefore, the expansion will converge well in the high temperature region. 

Other observables can also be calculated similarly. The quark number density 
nj is given by 

"/ 1 dlnZ d(p/T*) 

T3 VT*d{nf/T) d(n f /TY 1 J 

where T is fixed in the differentiations. When we define the light quark number 
density as n q = n u + n d , the susceptibilities of the light quark number density (x q ) 
and the strange quark number density (% s ) are given by 

X q _ ( d d \ n u + n d Xs _ d(n s /T 3 ) 

^ at.. lrr\ TlS ' rrO Q/.. /rp\ > I 4 ' 'J 



(4-8) 



T2 \d(f, u /T) d^a/T)) T3 ' d{ji a /T) 
respectively. The susceptibility of the isospin number can also be given as 
Xi f d d \n u -n d 



T 2 \d(fi u /T) d(fi d /T)J T3 



These quantities are expanded around jl = in terms of Cijk defined in (4-5). 
The trace anomaly is given by Eq. (2-3). The entropy density is given by s = 

T _1 + p — ^2jr fJ-fTif^j . The chiral condensate is defined by the derivative of InZ 

with respect to the quark mass. Taylor expansion of them can also be derived. 

4.2. Reweighting method and sign problem 

Another popular approach to finite density QCD is the reweighting method, 40 )~ 43 ) 
adopting the reweighting technique 44 ) ,45 ) to finite density QCD. Using the configu- 
rations generated at j2 = 0, expectation values at finite jl are computed by correcting 
the Boltzmann weight with the "reweighting factor" : 



Oxil f [detMf(» f )/detMf(0)] 

{0)<BRO) = v • (4'9) 

(n/tdetJlf/CM/VdetM/CO)]) 

The denominator is the ratio of the partition functions at finite jl and jl = 0, 

Z(j3,K,jl) _ / u detMffa f ) \ 

Z(/3,k,0) \11 detM/(0) / { ' 

Here, because det M* (/j,f) is complex, the reweighting factor 

n det Mf(n f ) 
detM/(0) 
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has a complex phase e td . When the fluctuation of 9 is larger than 0{ir/2), a reliable 
calculation of the numerator and denominator in Eq. (4-9) is difficult. This difficulty- 
is called the "sign problem (complex phase problem)". We actually encounter large 
fluctuations of 9 at large [if and/or large lattice volume. 

It is worth rewriting the denominator of Eq. (4-9) in terms of the distribu- 
tion function for P = —Sg/(6N S i te j3) (which is the plaquette if S g is the standard 
gauge action) and the logarithm of the absolute value of the reweighting factor 



F = In 



Uf (detAf'(/i/)/det Af'(O)) 



det M* (n f ) \ 



^ detM/(0) y/ 



(/3,k,0) 



=1 



(-) e F w (P,F;{3,K,0)dPdF, (4-11) 

P.F 



where wq is the distribution function for the phase-quenched system, 



w (P,F;P,FZ,jl) = J VU5(P - P[U])5(F - F[U]) 



JJdet M f (n f ) 
f 



a 6N site i3P 



(4-12) 



and {e %6 )p^F is the expectation value of the operator e ld at fl = with fixing the 
values of P and F to P and F: 



8{P - P) 5{F - F) e l 



p,f 



5{P - P) 5(F - F) 



(4-13) 



(/3,k,0) 



By measuring the histogram of P and F in a phase-quenched simulation, we 
can determine wq around the peak of the histogram. However, in the calculation 
of Eq. (4-11), precise information of wo is required around the maximum of the in- 
tegrant, which may deviate from the peak of wq due to the factor (e ld )p^e F . To 
calculate w§ in a wider range of (P, F) , we combine results of phase-quenched simu- 
lations at different points in the coupling parameter space, adopting the reweighting 
formulae for the phase-quenched theory. 13 - ) Further demonstration of such calcula- 
tion will be given in Sect. 6. 

For Eq. (4-11), we also need to estimate (e ld )p : F- Because the total distribution 
function is real and positive in finite-density QCD due to the charge conjugation 
symmetry, the imaginary part of (e ie }p i F must be averaged out when the statistics is 
infinite. Since the imaginary part is the source of the sign problem, we may remedy 
or mitigate the problem if we could find a method in which the imaginary part is 
removed and the real part is reliably estimated. In the next subsection, we show 
that it is useful to consider a cumulant expansion of {e ld )p^ in which ln(e )p,f is 
separated into real and imaginary parts. 13 )> 46 ) 

4.3. Cumulant expansion and Gaussian approximation 

For simplicity, let us consider the case of Nf degenerate quarks in the followings. 
A crucial step in handling the fluctuations in the phase 9 is to introduce an appro- 
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Fig. 6. The histogram of 6 in Nt = 2 QCD with improved Wilson quarks at (m, /m p ,T /T pc ) 
(0.65,0.94) (left) and (0.65,1.32) (right). 13 ' 



priate definition of 9 removing the ambiguity of complex phase with modulo 2tt. We 
uniquely define the complex phase by the Taylor expansion as 



0(/x) = iV f Im[lndetM(/x)] 



def. 



n=0 



(2n + 1) 



-Im 



9 2T,+1 (lndetM(/i)) 



9(/i/T) 



2n+l 



At=0 



(4-14) 



where the derivatives of IndetM can be unambiguously expressed in terms of M _1 
and derivatives of M. Note that the expectation value of 6 defined by Eq. (4-14) is 
not restricted to the range (— vr, tt), and the maximum value of \8\ is proportional to 
the volume of the system. The conventional complex phase in the range (— vr,7r) is 
recovered by taking the principal value of 6 with modulo 2tt. 

Typical results for the histogram of 9 are shown in Fig. 6 for JVf = 2 QCD 
with improved Wilson quarks, 13 ) where the power series in (4-14) is evaluated up 
to (/i/T) 3 discarding O ((/u/T) 5 ) terms*). We see that the width of the distribution 
becomes wider as /i/T is increased, indicating a severer sign problem. 

To mitigate the sign problem, we evaluate (e td ) p^f by the cumulant expansion: 46 ) 



p,F 



exp 



) c i(9 3 ) c , <<? 4 > c , i(9 5 ) c (0 6 ) c 



3! 



+ 



4! 



5! 



6! 



+ 



where {9 n 
/a6\ — It 



(4-15) 

c is the n th order cumulant: {9 2 ) c = (9 2 ) P , F , (9 4 ) c = (9 4 ) P)F - 3(9 2 ) PF , 
,„ lc = (9 6 ) P)F - 15(^ 4 )p 5 f(0 2 )p,f + 30(6> 2 )J, F , • • • . A key observation in handling 
the cumulant expansion is that (9 n ) c = for any odd n due to the symmetry under 
9 — > —9. This implies that, provided that the cumulant expansion converges, (e l9 )p : F 



*' In two-flavor QCD with p4 improved staggered quarks, 46 ' influence of the next order term 
was shown to be small at /x/T<2.5. See, e.g., Fig. 9 of Ref. 46). 
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is guaranteed to be real and positive and the sign problem is resolved. The sign 
problem is transformed into the convergence problem of the cumulant expansion in 
this approach. 

As shown in Fig. 6, we find that the distribution of 9 is well described by a 
Gaussian function up to fx/T ~ 0(1). The Gaussian distribution of 9 is observed 
also with improved staggered quarks in two-flavor QCD, 46 ) and was discussed in 
Ref. 47) too. The Gaussian distribution means that the cumulant expansion (4-15) 
is dominated by the lowest non-trivial order of (u), and thus the expansion is well 
converged. Corrections to the Gaussian distribution can be incorporated by taking 
higher order terms in the cumulant expansion. 

Here, we note that 9 = O(fifT) from Eq. (4-14). Therefore, if we take into 
account the cumulants up to the n th order, the truncation error does not affect 
the Taylor expansion up to O ((/i/T) n ). This means that the convergence of the 
cumulant expansion is closely related to the convergence of the Taylor expansion in 
\x. The Gaussian approximation is valid at least at small fi and the higher order 
cumulants may become visible at large fj,. The applicability range of the Gaussian 
approximation in jj, has to be checked for each cases by calculating higher order 
cumulants. 

We now argue that the range of applicability does not change with the system 
volume on sufficiently large lattices when the correlation length of the system is 
finite: 13 ) The expansion coefficients for 9 in Eq. (4T4) are given by traces of products 
of M _1 's and d n M/d(fj,/T) n, s over the spatial positions. For example, the first 
coefficient is given by the trace of Nf[M~ 1 (dM/d(fi/T))]. We note that the diagonal 
elements of this matrix are the local quark-number density operators [~ ■tp{x)'yo'ip{x)~\ 
at n = 0. When the correlation length of the local operators is finite and is much 
shorter than the system size, we may decompose the trace into a summation of 
independent contributions from spatially separated regions. The same discussion is 
applicable to higher order coefficients too. Then, the phase 9 may be written as a 
sum of local contributions from spatially separated regions, which are statistically 
independent with each other, i.e. 9 ~ ^2 x &x, where 9 X is the contribution from a 
spatial region labeled by x. The average of exp(z#) is thus 

(e»)»n(e*)=«p(EEs(^)' (4 ' 16 » 

x \ x n / 

This suggests that all cumulants (9 n ) c « J2 x (®x)c increase in proportion to the vol- 
ume as the volume increases, that is, (9 n ) c oc volume for any n, in contrast with 
a naive expectation that 9 n may be proportional to (volume)™ since 9 is propor- 
tional to the volume.*^ Therefore, while the width of the distribution, i.e. the phase 
fluctuation, increases in proportion to the volume, the ratios of the cumulants are 

This property of the higher order cumulants can be understood also from the effective poten- 
tial Veff(P, F) — — In wo (J 3 , F) — ln(e !e )(p F ), which will be studied in Sect. 6. Because V e a and In wo 
are proportional to the system volume, each term in the expansion of ln(e l9 ) (p,f) = E n ('°/ n ')^™)c 
should not increase faster than (volume) for any n. Otherwise, V e g becomes singular in the ther- 
modynamic limit. 
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independent of the volume — the higher order terms in the cumulant expansion are 
well under control in the large volume limit. 

Thus, the range where the cumulant expansion is applicable is independent of 
the volume. This is a good news in attacking the sign problem, which is known to 
become severer with increasing the lattice volume. Furthermore, we find that, when 
the system size is sufficiently larger than the correlation length, the distribution 
of ^/volume tends to a Gaussian distribution, since the distribution function of an 
average over many independent variables 9 X tends to a Gaussian function by the 
central limit theorem. 

4.4. Other approaches 

Besides the Taylor expansion method and the reweighting method, as well as 
combinations of these two methods, various methods have been proposed as alterna- 
tive approaches to study QCD at finite densities. An approach is to perform analytic 
continuation from simulations at imaginary chemical potentials. 48 ^' 49 - ) For a com- 
plex [if, Eq. (4-3) is generalized to \M* {Hf)y = jsMf(— [1*^)75. Therefore, when /j,f 
is purely imaginary, the Boltzmann weight is real, and we can simulate the system 
without the sign problem. Using results of simulations performed at imaginary /x's, 
information at small real [i can be obtained by an analytic continuation. The ana- 
lytic continuation is usually based on a Taylor expansion in terms of fx around fx = 0, 
i.e. we fit observables obtained at imaginary /x's with the Taylor expansion ansatz 
and extrapolate the resulting fitting function to a small real \i. Improvements of the 
analytic continuation procedure have also been discussed in Refs 50), 51) to obtain 
results in a wider range of real [i. Systematic high precision simulations in a wide 
range of the imaginary (i are required for a reliable analytic continuation. 

Another approach is to construct the canonical partition function Zq (T, N) by 
fixing the total quark number A or the quark number density. 41 )~ 43 )' 52 )~ 54 ) Using the 
canonical partition function, we can also discuss the effective potential as a function 
of the quark number. The relation between the grand canonical partition function 
Z(T,/j,) and the canonical partition function Zc(T,N) is given by 

Z(T,n)= (vU[&etM(n)] Ni e" 5 « Z c {T,N)e N ^ T (4-17) 

for the degenerate Af-flavor case. Because this is a Laplace transformation, Zq is 
obtained from Z by an inverse Laplace transformation. To find A that gives the 
largest contribution to Z, it is worth introducing an effective potential as a function 
of A, 

V cS (A, T,n) = — In Z C (T, N) - N ± = - N ± (4-18) 

where / is the Helmholtz free energy. V e g is useful to study the nature of phase 
transitions. 54 ^ 
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§5. Two flavors of improved Wilson quarks at finite density 

Adopting the Taylor expansion method and the Gaussian method discussed in 
the previous section, we made the first study of finite-density QCD with Wilson-type 
quarks. 13 ) We study two-flavor QCD with RG-improved gauge action and the clover- 
improved Wilson quark action. A systematic study of finite-temperature QCD with 
this action has been done by the CP-PACS Collaboration. 7 )' 9 ) The phase diagram 
at (i = as well as LCP's determined with m n /m p are given in Refs. 7), 9), 13). We 
extend the study to [i ^ 0. 

5.1. Taylor expansion for EOS 

We study the pressure and quark number densities defined by Eqs. (4-4) and 
(4-6) for the case of two- flavor QCD. Defining the quark chemical potential [x q = 
(fi u + Hd)/^ and the isospin chemical potential /i/ = (/j, u — fid)/2, the quark number 
and isospin susceptibilities are given by 

Xq _ d 2 (p/T±) xi _ d\p/T±) 

T 2 d{n q /Tf T 2 dim/T) 2 ' 1 } 

which measure the fluctuations in baryon and isospin numbers in the medium. 55 ) 
We study the isosymmetric case \i u = fi^ = /x, m = 0. The pressure is given by 

p ^ //ixn 1 AT3 d n lllZ 



J>(T)(!)", c n (T) = 



T 4 ^ ny >\TJ ' " v ' n\ Ni d(ti/T)> 



(5-2) 



Here, co(T) is the pressure at jl = and has been computed by the CP-PACS 
Collaboration. 7 )' 9 ) The susceptibilities are then expanded as 

^-Sc. + Mc. (£)' + ..., Se>) =2c , + 12cJ ^) 2 + ..., (5 .3, 

where 

j _ 1 AT t 3 cPlnZ(T, / u + / u / , / u-^ / ) 



n n\ Ni d(m/T) 2 d(n/T)n~ 2 - =0 • 1 " j 

In this study, we compute the Taylor expansion coefficients for the second and 
fourth derivatives in (5-2). This enables us to compute Xq an d Xi to the lowest 
non-trivial order in /i. 

5.2. Random noise method 

To evaluate the Taylor coefficients (5-2) and (5-4), we calculate 

d n IndetM 



V n = N i 



(5-5) 
fi=o 



up to n = 4. 13 ) We thus study 

_! dM \ 



V 2 = Ni 



f*, l d M \ f*, l dM », l dM \ 



(5-6) 
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etc., where 



x,y 



— K 



jl=0 



(1 - 14)U XA S X+Iy + (-lf(l + li)ul_ lA 5 x _ h 



(5-7) 



In terms of V n , we have c 2 = (N t /2N*){{V 2 ) + <Pf>}, = (N t /2N*) (V 2 ), 
Note that P n is real for even n and purely imaginary for odd n. 37 ) 

These traces, say trA, can be evaluated by "the random noise method". In 
this method, instead of calculating the diagonal elements An individually for all i, 
we calculate r^A-q for several random noise vectors rf. The contributions of the off- 
diagonal elements Aij {i / j) in this quantities are removed by averaging over the 
random noises, r)*rjj = 5ij. 

As is clear from the construction, the random noise method is effective when 
the contaminations from off-diagonal elements are small. Because the propagator 
(M~ 1 ) Xjy decreases rapidly with increasing \x — y\, the random noise method will 
work well to suppress small contaminations of spatially off-diagonal elements. On 
the other hand, the off-diagonal elements in the color and spinor indices are from 
the same spatial point and thus are not suppressed by \x — y\ — they are suppressed 
only by l/\/iV no ; sc , where N no i se is the number of noise vectors. This motivates us to 
apply the random noise method for the spatial coordinates only. For the trace over 
the color and spinor indices, we just repeat the calculation generating the random 
noise vectors for each color-spin index*). 

For a product of traces, the random noise vectors for each trace must be in- 
dependent. We compute such product by subtracting the contribution of the same 
noise vector from the naive product of two noise averages for each trace. This effec- 
tively increases the number of noises to N no - lsc • (iV no i sc — 1) for the products and thus 
suppresses their errors due to the noise method. 

In practice, N no { se can be small when the error due to the smallness of N no i sc is 
smaller than the statistical errors from the averaging over the configurations. The 
required number of noise vectors depends on each operator. In Ref. 13), we have 
tested the random noise method in the evaluation of T> n (with n = 1-4) . We find that 
T>\ has larger fluctuations under a variation of the noise vectors than T>2 etc., and 
the error in V\ dominates in the errors of C4 and c\ when we adopt the same AT noise . 
From this test, we choose AW = 100-400 for tr[(d" M /d{^ia) n )M- 1 ] (n = 1-4), 
and Ar no i se = 10 for other operators. 



*' Because a staggered-type quark does not have the spinor index at a spatial point, the number 
of off-diagonal elements is only 6 in the color 3x3 matrix, the contamination of off-diagonal elements 
is less serious. This is a reason that the random noise method is adopted more naively with the 
staggered-type quarks. However, with Wilson-type quarks, the number of the off-diagonal elements 
with similar magnitude in the quark matrix is f f times larger than the diagonal one. Therefore, the 
color-spinor index should be treated more carefully with Wilson-type quarks. 
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Fig. 7. Results of two-flavor QCD at finite density by the standard Taylor expansion method. 13 ' 
Results are obtained at m^/m p = 0.65. The truncation error of the Taylor expansion is 0(fi 6 ). 
To is the pseudo-critical temperature at [i — 0. Left: p-dependent contribution to the pressure, 
Ap/T 4 =p(^)/T 4 -p(0)/T 4 . Right: Quark number susceptibility Xi- 



5.3. Quark number densities and susceptibilities at /j, ^ 

Simulations are done on a 16 3 x 4 lattice in the fixed N t approach in the range 
/3 = 1.50-2.40 (T/Tp C ps 0.8-3 or 4 on LCP's corresponding to m^/nip = 0.65 and 
0.80, where T pc is the pseudo-critical temperature for each LCP. See Ref. 13) for 
details. 

Our results with the standard Taylor expansion method are shown in Fig. 7 for 
LCP at m^/rrip = 0.65. The left panel is for the finite density correction Ap/T 4 = 
p(fi)/T 4 — p(0)/T 4 of the pressure at finite \x. To is T pc at /i = 0. Recall that 
we have a crossover at T pc at /x = 0. The pressure changes more sharply as n is 
increased. When we increase //, Ap/T 4 becomes the same size as p/T 4 at n = 
around /x/T ~ 0(1). In the right panel of Fig. 7, we show the results of the quark 
number susceptibility x<j(aO at /i ^ 0. We see that a peak seems to be formed when 
we increase /i. However, in spite of various improvements in random noise estimators 
etc. as discussed in the previous subsection, the statistical errors due to the complex 
phase fluctuation of the quark determinant are still a bit too large to draw a definite 
conclusion about the peak. 

In order to suppress the errors from the phase fluctuation, we apply the cumulant 
expansion method discussed in Sect. 4.3. We compute the quark determinant ratio 
in Eq. (4-11) by the Taylor expansion up to 0(^i 4 ), 

N m ax -i 
n=l * ' 

-^V max -i 

0(/i) = iV f Im[lndetM( M )] « E (2n - i)] Imp 2n-i (pa) 2n -\ (5-9) 

n=l ^ 

with iV max = 2, and estimate the phase factor by the second order cumulant assuming 



F(n) = N f Re 



In 



det M(n) 
det M(0) 
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Fig. 8. The same as Fig. 7 but with the combined hybrid and Gaussian methods. 13 ' Left: Ap/T 4 
p(p)/T 4 — p(0)/T 4 . Right: Quark number susceptibility \q- 



10 - 



5 - 



Q 
0.5 





1 1 1 1 
H /r=i.2 




H /r=i.o 




n/r=o.8 




H /r=o.6 ~ 




n/r=o.4 




n/r=o.2 




n /r=o.o 






1 — 1 1 1 ■ 




i , i 



IT 1 



n /r=i.2 

n /r=i.o 

n/r=o.8 

H /r=o.6 

n /r=o.4 

H /r=o.2 

n/r=o.o 



r/r„ 




Fig. 9. Results of two-flavor QCD at finite density. 13 ' Left: Isospin susceptibility \i a t m^/m p 
0.65 by the standard Taylor expansion method. Right: Quark number susceptibility Xq 
m w /m p — 0.80 by the combined hybrid and Gaussian methods. 
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the Gaussian distribution of 6, {e ld ) ~ exp[— (6 2 )/2\. We moreover shift f3 from 
the simulation point /3n such that the statistical error due to the i^-integration is 
minimized in Eq. (4-11). We perform a fit of the resulting pressure (= free energy) 
in terms of \i to obtain Xq- The results for Ap/T 4 and Xq/T 2 are shown in Fig. 8. 
We find that the statistical errors are appreciably suppressed by these methods. 
We also note that, although the simulations at different T are independent, the In- 
dependences of Ap and Xq are smooth and natural. We can now clearly identify a 
sharp peak in Xq/T 2 that appears around T pc when (j,/T > 0(1). The peak becomes 
higher as fj, increases. These are consistent with the observations with staggered- type 
quarks and suggest a critical point at finite /i. 

In contrast with the peak in Xq, the isospin susceptibility shows no sharp peak, 
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as shown in the left panel of Fig. 9. This is in accordance with the expectation that 
Xi is analytic at the critical point since the iso-triplet mesons remain massive. 

Results at m^/mp = 0.80 are similar, but with milder peaks in Xq than those 
in Fig. 8, as shown in the right panel of Fig. 9. This may be explained in part by 
the expectation that the critical point locates at larger fx because the quark mass is 
larger than that for m^/nip = 0.65. See Ref. 13) for more discussions. 

§6. Histogram method and QCD phase structure at zero and finite 

densities 

In a study of the QCD phase structure shown in Figs. 1, and 2, identification 
of first-order transition region is quite important. Among several methods to study 
the nature of phase transitions, the probability distribution of physical observables 
provides us with the most intuitive way: The probability distribution function is de- 
fined as the generation rate of configurations with fixed expectation values of physical 
observables. Double or multiple peaks in the probability distribution function for 
observables which are sensitive to the phase, such as the energy density, chiral order 
parameter, etc., give a signal of first order phase transition. A problem is that, in 
order to trace the variation of the shape of a probability distribution function, we 
need a statistically reliable data of the distribution function in a wide range of the 
expectation values. In the identification of a first order transition, a correct evalua- 
tion of the double-peak distribution requires a quite long simulation with sufficiently 
many flip-flops among different phases [see, e.g., Ref. 56)]. This is computationally 
quite demanding with dynamical quarks. 

Here, we note that the calculation of the probability distribution function is 
required also for the reweighting method for the components of the action. 44 )' 45 ) 
Therefore, when we adopt observables which are the components of the action, the 
computation of the distribution function at different points in the coupling parameter 
space is straightforward by the reweighting method. 46 ) This helps us to obtain the 
probability distribution function in a wide range of expectation values. Therefore, 
the method is quite powerful to study the location of first order transitions. 

We apply the method to explore the phase structure of QCD. In this section, we 
introduce the method, which may be viewed as a variant of the histogram method or 
the density of state method, 57 ) and test it in the heavy-quark region of QCD. 14 )' 15 ) 
We then present our on-going project to study finite density QCD with light dynam- 
ical quarks by combining the histogram method with phase-quenched simulations. 16 ) 

6.1. Histogram method 

As a demonstration, let us consider the simplest lattice QCD: the combination 
of the plaquette gauge action with unimproved Wilson quarks. 



S — S g + S q , 
S g = -6N site f3P, 



(6-1) 
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N t r 3 

f=i x y M =i 

/=1 x,S/ 

where iV s ite = iV s 3 x iV t is the lattice volume and 

P = W~ E \ Re tr K^+A,^ +i >,X J (6-3) 

Sltc X, [1<V 

is the plaquette. Note that M does not depend on /?.*) 

Denoting the values of the operators (P, • • • ) as (P, • • • ), the probability distri- 
bution function for (P, • • • ) is defined by 

w{P,--- \p,H,p) = !vU5(P[U]-P} ■ ■ - JJdet M(n f , n f ) e 6Nsi ^ p 
J f 

= e 6N sitcf sp fT>U6(P[U]-Py--Y[debM(Kf,Hf), (6-4) 
J f 

where "• ■ • " in the r.h.s. means the product of delta functions for other operators. 
We now define the effective potential as 

V cS (P,--- ;P,K,p) = -hiw(P,--- \p,K,fi. (6-5) 

Note that P represents the freedom of gauge internal energy, and thus should be 
sensitive to the phase structure of the system. 

A useful property of the plaquette distribution function and the effective poten- 
tial is 

w(P, ■■■ ;P,K,p) = w(P,--- ; /3o, k, U) eW-M N «** p , (6-6) 
V cS (P, ■ ■ • ; P,k, fx) = V cS (P, ■ ■ ■ ; A), k, fl) - 6(/3 - f3 )N sitc P. (6-7) 

We thus find that 

^r(P,--- ;/3,«,m) = ^(P,-" ;/3o,K,/I)-6(/3-/3o)iVsitc, (6-8) 
and d 2 V e ff/dP 2 is independent of /3. 



*^ When we consider improved gauge actions such as (3-2), we replace P with the operator 
appearing in the gauge action: P =' —S g /(6N S it e f!l). On the other hand, when M depends on /3 
as in the case of improved quark actions, more careful treatments are required. See discussions in 
Sect. 6.5. 
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The k and /7-dependences of V e e can also be computed by the reweighting method 
as follows: 

V e s(P,--- I = V eS (P, ■■■ ;P,Ko,th) - ]nR(P,--- ; Ft, fi, kq, fio), (6-9) 

where the reweighting factor R is evaluated as 

R{P, -- ;K,fj,,Ko,fj, ) - 



w(P,--- ;P,K ,fo) 

JVU5(P -P) -Uf det M ( K f, M/) e 6/37Vsit ^ 
JVU5(P - P) ■ ■ ■ H f det M (kJ, /xj) e 6 ^^ 
/P^P-^-.-ri/detM^/x/) 
/ PC/5(P - P) ■ ■ • 11/ det M («/> M/) 
a^p p^i rr dctM ( K f^f) 

oK^-n- ■■ii/ditM(ra) // 

- 1 - 1 /(K(),^o) 

«5(P-P)...) ' 
/ n detM(K^\ 

V 7 7 J ' (Xo,ik);P,- 

Note that P is independent of (3 and thus can be evaluated at any f3. By adjusting 
and combining (3, we can obtain precise values of R in a wide range of P, • • • . 

6.2. (JCP in the heavy-quark region 

We first test the method in the heavy-quark region. As shown in Fig. 2, we have 
the first order deconfinement transition of the SU(3) Yang-Mills theory in the limit 
of infinite quark masses. The transition is expected to turn into a crossover when 
we decrease the quark masses. We study the boundary of the first order region by 
the histogram method. To take advantage of light computational costs in quenched 
simulations, we choose = /j/j = (/ = 1, ■ ■ ■ , Nf). 

For simplicity, let us consider the case of degenerate quarks: kj = k, [if = n 
(/ = 1, • • • , Nf). Generalization to non-degenerate cases is easy. 

To the lowest order of the hopping parameter expansion, the quark determinant 
ratio appearing in the r.h.s. of (6-10) is evaluated as 



exp 288N sitc K 4 P + W*2 Ni+2 K Nt {cosh f2 R + isinh (£) A} 



det M(k,(i) 

detM(0,0) 

(6-11) 

where 

s x 

is the Polyakov loop, and J7r = Re.!? and Q\ = Imi? are its real and imaginary parts. 
We note that the contribution of P in the r.h.s. can be absorbed by a shift of the 
gauge coupling /3 — > (3 + 4&N{K 4 . 
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Fig. 10. Effective potential at ji — as a function of P in the heavy-quark region. 14 ' Left: 
dV a a(P; /3,0) jdP in the heavy quark limit. Data obtained at five different values of /3 in the 
range 5.68-5.70 are converted to /3 = 5.69 by (6-8). Right: V e a(P; j3, k) for k > 0, where /3 is 
adjusted such that the two minima of V e g have the same depth and the constant term of V c s is 
adjusted such that V e « = at the central peak point P poa k- 



The term proportional to Q\ induces the complex phase at // ^ 0, which is the 
origin of the sign problem at large \x. 

6.3. Results at [i = in the heavy-quark region 

We are now ready to calculate the effective potential. Let us first study the case 
H = 0. 14 )' 15 ) fn this case, the complex phase term is absent in (6-11). 

The double- well nature of the effective potential is clearly seen when we consider 
the effective potential for P only: 

V cS (P; (3, k) d = - In J w(P, f2 R ; [3, k) dQ K . 

Our results for dV e ^{P] (3, 0)/dP at k = are shown in the left panel of Fig. 10. 
Using (6-8), we shift the results obtained at five different /3's to /3 = 5.69. With 
different (3, the range of P in which we can reliably obtain V e s is different. We find 
that the results of dV e g (P; /3, 0)/dP at different /3's are smoothly connected with 
each other by (6-8). We can thus obtain accurate values of dV c fi (P; f3, 0)/dP in a 
wide range of P. 

Similarly, we calculate dV e g(P; /3, K)/dP at k > by using the reweighting factor 
R and (6T1). We then integrate dV e s {P; f3, K)jdP in P to get V e ff(P; f3, k) as shown 
in the right panel of Fig. 10. We find that the double- well structure of V e ff(P) 
becomes weaker with increasing k, and eventually disappears at finite k, say K cp . 
Examining the shape of V e g more closely, wc obtain K C p — 0.0658(3) (I 4 ! ) for two- 
flavor QCD in the lowest order of the hopping parameter expansion on an Nt = 4 
lattice. 14 ) 

The argument can be easily generalized to the case of non-degenerate quark 
masses. Our results for the critical point K cp in 2 + 1 flavor QCD is shown in 
Fig. II. 14 ) The top-right corner of Fig. 2 corresponds to this plot rotated by 180°. 
Our results are consistent with that obtained in an effective model 60 ) and with a 
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Fig. 11. The phase boundary separating the first order transition region and crossover region in 
the (reud, K B ) plane at fi — in the heavy-quark region. 14 ' 

recent study using the hopping parameter expansion. 61 ) 

The expression (6-11) suggests us to adopt Q R as an additional operator for the 
effective potential V e Q. Then, the f2 R in the r.h.s. of (6-11) is simply replaced by its 
expectation value and the reweighting factor R is just a given function of P and 
1?r in this case. We have 



° 1 ' ;ff (P, Q R ; fi, k) = ^(P, f2 R ; fi ) - 67V site (fi + 48N { k a - fi ) 



dP 

dV eS 



dP 

^(P,f2 R )-3N^ +2 N iK N \ 



(6-13) 
(6-14) 



df? R df2 R 

where Vo is the effecvtive potential in the heavy quark limit (SU(3) pure gauge 
theory). The argument /3 in dV e ft/df2 R and dVo/df2 R is omitted in (6-14) since they 
are independent of due to (6-7). Note that, besides known overall constants [the 
last terms in (6-13) and (644)], the dependences of dV e d/dP and dV e d/df2 R on P 
and J?r are independent of fi and k. 

Because f2 R represents the freedom of heavy-quark free energy, we expect it be 
sensitive to the phase structure of the system. Around the first order transition point, 
we will have a double- well structure of V e g in the two-dimensional plane (P, f2 R ). To 
study the phase structure, it is useful to examine the curves dV e g/dP = and 
9V e f[/df2 R = 0. From (6-13) and (6T4), these curves at different (fi, k) corresponds 
to different contour curves of dVo/dP and dVo/df2 R . A contour plot for dVo/dP 
and 8Vo/df2 R is given in Fig. 12. When the curves dV e s/dP = and dV e g /df2 R = 
cross at only one point, we have just one minimum of V e g. In this case, we 
have no first-order transitions around this (fi,K). On the other hand, when we 
have three intersection points, we have two minima and one saddle point, implying 
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Fig. 12. Contour plot of dVo/dP (blue curves) and 9Vo/9J7r (red curves) at /i = in the heavy- 
quark region. 15 ' Values of [3* = /3 + 48N{K 4 and k for the corresponding curves of dV e s/dP = 
and dVctt/ dfin. = are also given. 



the existence of the first order transition. In particular, from the merger of three 
intersection points, we can determine the critical point where the first order transition 
line terminates. From Fig. 12, we find that the S-shaped curve of dV c s/dQYi = 
leads to three intersection points at small k, and that the S-shape becomes weaker 
with increasing k. A preliminary estimate for the critical point is K cp ~ 0. 0690(7), 15 -* 
shown by thick contour curves in Fig. 12. We are currently testing a refinement of 
the method to extract smoother contour curves. 

To examine the quality of the hopping parameter expansion, we have studied 
the effect of the next-leading order terms in the evaluation of the critical point 
K cp . 58 )' 59 ) To this order, we need to incorporate K 6 -loops with length six and gener- 
alized Polyakov loops with length Nt + 2 in the quark determinant ratio. The effects 
of K 6 -loops may be absorbed by a shift of (3. Examining the effects of generalized 
Polyakov loops, we find that At cp shifts only by about 3% on the Nt = 4 lattice due to 
the next-leading order terms. We also find that the contribution of the next-leading 
order terms becomes comparable to that of the leading order terms at k ~ 0.18. 
Because K cp is much smaller than this, we conclude that the hopping parameter ex- 
pansion is well valid up to k ~ K cp . Accordingly, an estimation of the pseudo-scalar 
meson mass around K cp leads to m n w T/0.023 for Nf = 2, 14 ) i.e., m n ~ 7-9 GeV 
with T 



160-200 MeV. Thus, K cp locates well in the heavy-quark region. 



6.4. Results at \i ^ in the heavy-quark region 

In order to extend the study to finite densities, we have to calculate the reweight- 
ing factor due to the complex phase, 




with = 3N*2 Nt+2 N f \f2i, A = K Nt sinh (/i/T) . (6-15) 
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where (• • • )p,n n is the expectation value at fixed P and 17r in quenched QCD. When 

9 fluctuates a lot at large fi, a reliable estimation of (e )p,i? R becomes difficult (the 
sign problem).*^. 

Before evaluating (e t0 )p j Q Ii , let us consider the case of phase-quenched finite 
density QCD, in which the complex phase term is removed in the quark determinant. 
In two-flavor QCD, this corresponds to the case of the isospin chemical potential, 
= — /Xd = fJ>i- From (6-11), we find that, after shifting (3 — > j3 + 48iVf k 4 , the effects 
of fij are just to further shift k — > k cosh 1 ^' (fij /T). Therefore, we have 

«cp(Mi) = MOVcosh 1 /^//^) (6-16) 

for the critical point in the phase-quenched QCD to the lowest order of the hopping 
parameter expansion, where n cp (0) is the critical point at fi u = fid = 0. Note that, 
with increasing the critical point approaches towards k = where the hopping 
parameter expansion is reliable. 

We now compute the effect of the phase. To evaluate (e ie )p^ H , we adopt the 
Gaussian approximation with the cumulant expansion 13 )' 46 ) discussed in Sect. 4.3. 
In the heavy-quark region, we study {9 2n ) c = (3N*2 Nt+2 N { A) 2 " {f2? n ) c . We find 

that (0 2 ) c 3> (v )c around the critical point. 15 ) This confirms the validity of the 
Gaussian approximation. We thus have 



2\ 



^ A (f, + ^ - a) + (3^7^)° gig, ( , 17) 

_ H . uv~»* ( £) + sas^!^. ( , 18) 

When the last term in (6-18) modifies the S-shape of the curve dV e s / '(9J?r, = shown 
in Fig. 12, k cp (/j,) deviates from nf. p (fJ-)- By evaluating these quantities, we however 



*' If we could treat Qi as a variable for V c s too, the reweighting factor (e lS ) is just a given 
function of 0% (or equivalently 8). However, when the fluctuation in 8 is large, then it is difficult 
to determine a reliable V e g unless quite a high statistics is accumulated. This is a rephrasing of the 
sign problem. 
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find that the contribution of the last term in (6-18) is at most about 3% of the second 
term around K^. p (fi) even in the large (i limit. Therefore, K cp (fi) is indistinguishable 
from «c p (/x) within the current statistical errors: 15 ) 

Mm) ~ M°)/ coshl/JVt Wr) (6-19) 

Generalization of this result to non-degenerate cases such as the 2 + 1 flavor QCD 
(fi u = fid = Hud) is straightforward: 

2 Kd(fI)} Nt cosh(fi ud /T) + [n,(j£)] N * cosh^/T) « 2 [k% =2 (0)] n *. (6-20) 

Critical surfaces for the symmetric case fi u = fid = fi s = I 1 and a more realistic case 
of fi u = fid = fJ-ud and fi s = 0, that may be realized in heavy ion collisions, are shown 
in Fig. 13. 

6.5. Phase- quenched simulations towards the physical point 

We are challenging to apply the histogram method to explore the phase structure 
of finite-density QCD in the light quark region. For the sake of notational simplicity, 
we consider the degenerate case in this section too. Generalization to non-degenerate 
cases is straightforward. 

When quarks are light, the Polyakov loop Q plays no more a decisive role in the 
dynamics of the system. We thus consider detM itself as the additional operator 
for the effective potential, det M represents the freedom of quark internal energy, 
and thus should be sensitive to the phase structure of the system. We denote the 
absolute value and the complex phase of the quark determinant as follows: 

F(f3,K,n) = N{ \n\det M (p, K,fi)/ det M(f3,K,Q)\ 

~d In det M ((3, n,fi' 



N t / M Rc 
Jo 



dfi' 

6(P, k, fi) = N f Im [In det M((3, k, fi)\ 

~d In det M(l3,K,fi r ) 



Jo 



dfi', (6-21) 



dfi', (6-22) 



dfi' 

where d(lndetM)/<V = t^M^^M / dfi')] can be evaluated by the random noise 
method discussed in Sect. 5.2. Note that 6 is uniquely defined in the range (— oo, oo), 
as discussed in Sect. 4.3. The distribution function and the effective potential for P 
and F are now given by 

w(P,F;i3,K,fi) = JvU5(P - P)S(F - F) e i§ \det M(f3, k, fi)\ Nc e ^ N ^ p 

w (P,F;(3,K,n), (6-23) 



(0:P,k,h);P,F 

V eS (P,F;(3,K, f i) = V (P,F;/3,K,fi)-\n(e i < j ) (6-24) 

\ / (0:p,K,n);P,F 

where wo and Vq = — In wo are the distribution function and the effective potential for 

the phase-quenched system, and the expectation value (e i6 ')(o: ( g iKiAt ) ; p i _F is evaluated 
with fixed P and F in the phase-quenched simulation. 
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Fig. 14. Results of phase-quenched simulations in two-flavor QCD with RG-improved gauge action 
at f) = 1.5. 16 ' Left The effective potential Vo(F) at fj,/T — 2.0 evaluated at three different 
simulation points. Right: The distribution of the phase of the quark determinant at fi/T = 2.4. 
The dashed curves are the fitted results with the Gaussian function. B% is the fourth-order 
Binder cumulant normalized such that B% = 3 for the Gaussian function. 



In the following, let us consider a simpler case where the quark matrix M can 
be treated as independent of /3. E.g., in a study of the phase structure at a value of 
/3, we may treat f3 in M as fixed to that value. Physical properties such as the phase 
structure will not be affected by this procedure*). Then, (e^W^^vp^ becomes 
independent of j3 and can be evaluated at any /3 to cover a wider range of P and F. 

Because the phase-quenched simulations in two-flavor QCD correspond to the 
case of isospin chemical potentials, a comment is in order about the influence of 
the pion-condensed phase which exist at large isospin chemical potentials. &2 > In the 
pion-condensed phase, (e ld ) is expected to vanish by model studies. 63 )> 64 ) According 
to (6-23), this means that the configurations in the pion-condensed phase have no 
contributions to the physics of phase-unquenched QCD — w and V c s are dominated 
by phase-quenched configurations out of the pion-condensed phase, and we only need 
to generate configurations outside the pion-condensed phase. 

We test the method in two-flavor QCD adopting the RG-improved Iwasaki gauge 
action (3-2) and the clover-improved Wilson quark action (3-3). According to the 
footnote in Sect. 6.1, P is identified as 

Sltc X Kfl>U fl,U ) 

in the present case. We perform simulations at m w /m p ~ 0.8 on an 8 3 x 4 lattice with 
a non-perturbatively estimated csw- in the present study of the phase structure at 
each {}, we treat csw a s a constant independent of j3. 

In the left panel of Fig. 14, we show our results of Vq as a function of F (after 

*' When we calculate observables as functions of /3, we need to incorporate the effects from the 
/3-dependence in M. See 65) as a trial to incorporate the /3 dependence of csw in the reweighting 
factor. 
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integrating out P). The results of simulations at three values of /jlq/T are translated 
to jji/T = 2.0 using a reweighting formula for Vq in the //-direction: 

V (P,F;P,K,fi) = Vo(P,F;P,K,fio)-]nRo(P,F;K,n,fi ) (6-26) 
def. w (P,F;/3,K,fi) 



R (P,F;k,h,ho) 



w (P,F;f3,K,fi ) 

jVU5(P - P)S(F - F) | detM(n,(i)\ Nt e 6(3N ^ p 
fVU5(P - P)5(F - F) | det M(k,ho)\ n ' e*t )N *** p 
fVU5(P - P)S{F - F) | det M(k, fi)\ N f 
f VU5(P - P)S{F - F) | det M(k, n )\ Nt 
det M(k,h) Nf * 



det M{k,hq) 



(6-27) 

(0:^,^0);^,^ 



where Rq does not dependent on /3. From this figure, we confirm that the data 
obtained at different simulation points form a smooth Vo- We can thus obtain precise 
Vq in a wide range of F. 

To calculate (e ld ) (o-.k,h)-,p,f > we adopt the cumulant expansion method. Typical 
result for the distribution of 6 is shown in the right panel of Fig. 14. We find that 
the distribution is well approximated by a Gaussian function. We also note that the 
second-order cumulant increases with increasing fi, while the forth-order cumulant is 
consistent with zero within the statistical error, though the statistical error increases 
with ju. 16 ) Therefore, the cumulant expansion is well controlled by the leading term 
and we may reliably evaluate the complex phase factor even at these relatively high 
values of fi. A project towards clarification of the phase structure at the physical 
point is under way in this direction. 



§7. Heavy-quark free energy 



Finally, we study the heavy-quark free energies and screening masses in QGP. 
The free energies for static quark-antiquark and two static quarks characterize inter- 
quark interactions in QGP, and their Debye screening masses describe the thermal 
fluctuation of quarks and gluons in QGP. In a phenomenological model, they are 
relevant to the fate of heavy-quark bound states such as J/ij; and T in QGP created 
in the relativistic heavy-ion collisions at RHIC and LHG 66 )' 67 ) On the lattice, studies 
in quenched QCD 68 ) -71 ) and in full QCD with staggered-type quark actions 72 ) -75 -' or 
with Wilson-type quark actions 13 ^' 19 )~ 21 )> 76 ) have been made. Comparisons with 
analytic studies 77 )' 78 ) have also been attempted. 

Heavy-quark free energies on the lattice are defined through the correlation 
functions of the local Polyakov line operator J?(x) d =' IXr^i ^(x,t),4- Note that the 
trace over the color indices is not taken in i?(x). With a gauge- fixing, we define the 
free energy F R in various color channels i?: 79 ) 

e -F\r,T^)/T = I(Tr^t(xMy)), (7-1) 

3 
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Fig. 15. Heavy-quark free energies in two-flavor QCD for color singlet and octet QQ channels (left) 
and color anti-triplet and sextet QQ channels (right) obtained at m^/m p = 0.65 on a 16 3 x 4 
lattice. 19) The free energies are normalized such that they vanish at large distances. 



-FHr^/T = l (Tr ^t(x)1YJ7(y)) - l<Tr^( x)i 7(y)>, 



-F e (r,T,ii)/T 



-F 3 (r,T,»)/T 



l<Trtt(x)Trrt(y)> + l(Tri?(x)i7(y)), 
i(Tr^(x)Tri7(y)) - ^(TVfi(x)fi(y)), 



(7-2) 
(7-3) 
(7-4) 



where r = |x — y|, and R = 1 for the color singlet QQ channel, 8 for the color octet 
QQ channel, 3* for the color anti-triplet QQ channel, and 6 for the color sextet 
QQ channel, respectively. To preserve the free energy interpretation of F R by the 
transfer matrix theory, the gauge-fixing procedure should not include the temporal 
links. We thus adopt the Coulomb gauge. 

Above T pc , we also introduce normalized free energies V whose constant terms 
are adjusted such that they vanish at large distances r — > oo. This is equivalent to 
defining the free energies by dividing the r.h.s. of Eqs. (7-l)-(7-4) by (Trl7(x)) 2 . 

7.1. Heavy- quark free energies in two-flavor QCD 

We first compute the free energies (7-l)-(7-4) in two-flavor QCD at zero and 
finite densities. 19 ^ We consider the case fJ- u = = We use the gauge configu- 
rations generated for the studies discussed in Sect. 5. As mentioned in Sect. 5.3, 
the configurations were generated on a 16 3 x 4 lattice on LCP's corresponding to 
m n /m p = 0.65 and 0.80. We thus adopt the fixed- N t approach for this study. 

7.1.1. Heavy-quark free energies at /i = in two-flavor QCD 

The heavy-quark free energies at /i = are shown in Fig. 15 for m v /m p = 0.65 
and T > T pc . Results for m n /m p = 0.80 are similar. We find that the inter- 
quark interaction is "attractive" in the color-singlet and antitriplet channels and is 
"repulsive" in the color octet and sextet channels. We also see that, irrespective of 
the channels, the inter-quark interaction becomes rapidly weak at long distances as 
T increases, as expected from the Debye screening at high temperatures. 

We find that the heavy-quark free energies in the high temperature phase are 
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Fig. 16. The effective running coupling a e e(T) (left) and Debye screening mass m£>(T) (right) in 
two-flavor QCD, obtained at m w /m p = 0.65 on a 16 3 x 4 lattice. 19 ' 



well described by the screened Coulomb form, 

V R (r,T) = C R ^P-e- m ^ r , (7-5) 

where a e s(T) and mo{T) are the effective running coupling and Debye screening 
mass, respectively. From the fits, we note that the color-channel dependence of 
the free energies can be absorbed by the kinematical Casimir factor inspired by the 
lowest-order perturbation theory: 

Ci = ~o> ^8 = c' ^6 = ^3* = --, (7-6) 

This Casimir dominance at high temperatures was reported also in quenched studies 
using the Lorentz gauge. 69 ) With the Casimir factors, a e g(T) and itid(T) are well 
universal to all color channels, as shown in Fig. 16. The magnitude and the In- 
dependence of the Debye mass are also consistent with the next-to-leading-order 
thermal perturbation theory. 19 ) 

7.1.2. Heavy-quark free energies at fi ^ in two-flavor QCD 

The Taylor expansion of V R with respect to fJL q /T is given by 

V R (r,T, N )=v*(r,T)+v R (r,T) (^) +v R (r,T) [^.)\ ^\ (7-7) 

where concrete forms of the Taylor coefficients v R are given in Ref. 13). 

The color singlet and octet channels do not have the odd orders in the Taylor 
expansion since the QQ free energies are invariant under the charge conjugation. 
v R is the normalized free energy at \x = shown in Fig. 15. Results for the lowest 
non-trivial order are shown in Fig. 17. See Ref. 13) for higher orders as well as those 
at m-n/mp = 0.80. From these figures, we find that, both around T pc and at higher 
temperatures, the sign of v R is the same with that of v R , whereas the sign of a v R 
is the opposite of that of v R : 

v R -v R >0 (only for QQ free energies), v R -v R <0. (7-8) 
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Fig. 17. Taylor coefficients for heavy-quark free energies in two-flavor QCD, obtained at m^-/m p = 
0.65 on a 16 3 x 4 lattice. 13 ' Left: for color singlet and octet QQ channels. Right: V\ for 
color anti-triplet and sextet QQ channels. 



Because is absent for QQ free energies, this means that, in the leading-order of 
fj, q , the inter-quark interaction between Q and Q becomes weak at finite [A q , while 
that between Q and Q becomes strong. In other words, QQ (QQ) free energies are 
screened (anti-screened) by the internal quarks induced at finite \i q . 

Taylor expansion coefficients for a e ff(T, /j) and m£)(T, //) can be computed sim- 
ilarly. We find that the leading correction in mo (T, /i) due to finite fi is larger than 
a prediction of the leading-order thermal perturbation theory. 13 ) 

7.2. Heavy-quark free energies in 2 + 1 flavor QCD 

We now extend the study to the more realistic 2 + 1 flavor QCD. 21 )' 34 ) As 
discussed in Sect. 3, we adopt the fixed-scale approach for 2 + 1 flavor QCD. We thus 
vary T by varying Nt with fixed coupling parameters. We use the finite-temperature 
configurations generated in Sect. 3 to compute the heavy-quark free energies at 
m n /m p ~ 0.63 and m^ ss /ra^ ~ 0.74. In this study of 2 + 1 flavor QCD, we restrict 
ourselves to the case \i = 0. 

A good feature of the fixed-scale approach is that the renormalization factors, 
which are determined on a zero-temperature lattice depending on the coupling pa- 
rameters, are common to all T's in the fixed-scale approach, because the coupling 
parameters are fixed for all T's. This is so also for the heavy-quark free energies. 
Therefore, we can directly compare the bare free energies at different T's in the 
fixed-scale approach. 

Our results for the bare free energies are shown in the left panel of Fig. 18. 
Plotted are the data for the color singlet QQ channel at various T's in the high tem- 
perature phase, together with the zero-temperature static quark-antiquark potential 
V(r) defined by the Wilson loop operator: 

V{r) = - lim Ir 1 In (w[ 4 x£ )] . (7-9) 

We find that singlet free energy T 1 (r, T) at any T converges to V(r) at short dis- 
tances. This is in accordance with the expectation that the short distance physics is 
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Fig. 18. Free energies of static quarks in 2+1 flavor QCD at finite temperatures with the fixed-scale 
approach. 21 ' The scale was set by the Sommer scale ro = 0.5 fm. Left: Bare free energy in the 
color-singlet channel. The static quark-antiquark potential V(r) at T = has been calculated 
by the CP-PACS and JLQCD Collaborations. 28 ' The fit result of V(r) by the Coulomb + 
linear form is shown by the dashed gray curve. The arrows on the right side denote twice the 
single-quark free energy 2Fq. Right: Normalized free energies for color singlet and octet QQ 
channels. 



insensitive to temperature. With the fixed-scale approach, we directly confirm that 
the theoretical expectation is actually satisfied on the lattice. *) 

At large distances, i ?1 (r, T) departs from V(r) and eventually becomes flat due 
to Debye screening. On the right edge of the left panel of Fig. 18, we show twice the 
single-quark free energy defined by 2Fq = — 2Tln(Tri?(x)) at each T by the arrows 
with the same color. We confirm that i ?1 (r, T) converges to 2Fq(T) quite accurately. 

By subtracting 2Fq, we obtain normalized free energies shown in the right panel 
of Fig. 18 for the QQ channels. See 21) for the results in the QQ channels. Performing 
fits with the screened Coulomb form (7-5), we confirm the Casimir dominance (7-6) 
as in the case of two-flavor QCD. 

7.3. Gauge-independent screening masses 

The Debye screening masses and the effective couplings computed in the previous 
subsections are dependent on the choice of the gauge. Therefore, their physical 
meanings are not quite clear. In Ref. 20), we have proposed a gauge-independent 
definition of screening masses for electric and magnetic channels. 

Under the Euclidean time-reflection 1Z and the charge conjugation C, the gluon 
fields transform as, 

Mt,x) ^ ^(-r,x), A 4 (r,x) ^ -A 4 (-r,x), A„(r,x) ^ -A*(r,x). (7-10) 

We call an operator magnetic (electric) if it is even (odd) under 1Z. It is natu- 
ral to extract magnetic and electric properties by decomposing observables using 

*' In the case of the conventional fixed- Nt approach, different renormalization factor is required 
at each T. In early studies, the insensitivity at short distances was just assumed and used to adjust 
the constant term of F 1 (r, T) at different T. 70) ' 72) In more recent studies, the renormalization 
factors are computed by a series of zero-temperature simulations at the same coupling parameters 
as the finite temperature simulations. See, e.g., Ref. 80). 
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Fig. 19. Comparison 20 ' of the screening ratio, m E _/m M+ , with predictions in the dimensionally- 
reduced effective field theory (3D-EFT) 82 ' and TV = 4 supersymmetric Yang-Mills theory. 83 ' 

these symmetries: 81 ^ Under these transformations, the Polyakov line operator f?(x) 
transforms as 



ft(x) ^> tt f (x), l?(x) ^ f2*{x) 



(7-11) 



Magnetic (electric) Polyakov line operator can be defined as IZ-even (7£-odd) part of 

^ M (x) ee 1 (fl(x) + fl+(x)) , J7b(x) ee ~ (l?(x) - !?t( x )) , (7-12) 

which can be further decomposed into C-even and odd parts as 

4l±(x) = i (4l(x) ± i^(x)) , 4s±(x) = ~ (^e(x) ± J2 E (x)) , (7-13) 

where ± stands for even or odd under C. Note that Tr,i?M- = TV/2e+ = and 
Trj?M+ (Trf?E-) is nothing but the real (imaginary) part of Trfi. 

Then the magnetic (electric) screening mass is extracted from the long-distance 
behavior of generalized gauge- invariant Polyakov loop correlation functions, 



C M+ (r,T) = (Tr^M+(x)Tr^ M+ (y)) - (TW?) 
C E _(r,T) = (Tr^ E _(x)Tr^ E _(y)\ > B - 



->A- 



r— >oo 
-m E _(T)r 

rT 



-m M +{T)r 

rT 



(7-14) 



where r = |x— y | and (Tri?) is real due to the C symmetry. Note that the conventional 
gauge-invariant Polyakov loop correlation function is given by 



C n {r,T) = (Trl?t( x )Tr%)) - (Tri?) = C M+ (r,T) - C E _(r,T). (7-15) 



Using the configurations generated for 19), we computed these screening masses 
in two-flavor QCD in the high temperature phase. We find that (i) Cm+ and Cg_ 
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have opposite sign, and (ii) Cm+ has larger magnitude and longer range than Ce- 
at long distances. The latter implies that the the conventional Cq is dominated 
by the magnetic sector at long distances, and thus m n (T) ~ rriM+(T). We also 
find m,M+(T) < ms-iT). A comparison with a dimensionally-reduced effective field 
theory 82 ^ and an M = 4 supersymmetric Yang- Mills theory with AdS/CFT corre- 
spondence 83 ) lead to good agreements of m E _/m M+ for 1.5 < T/T pc < 3, as shown 
in Fig. 19. 20 ) Further study is needed to clarify the meanings and implications of 
these results. 

§8. Summary 

The WHOT-QCD Collaboration is pushing forward a series of projects to clar- 
ify the phase structure and thermodynamic properties of the QCD matter at finite 
temperatures and densities, mainly adopting improved Wilson quarks. Wilson-type 
quarks do not have the theoretical unclearness of the staggered-type quarks, but 
require more computational resources. Thus, development and application of vari- 
ous computational techniques are mandatory. We have developed the T-integration 
method to calculate the equation of state in the fixed-scale approach, the Gaus- 
sian method based on the cumulant expansion to tame the sign problem, and the 
histogram method combined with the reweighting technique to explore the phase 
structure of QCD. By adopting them, we have made a series of studies in QCD at 
finite temperatures and densities with improved Wilson quarks. In particular, we 
have carried out the first study of finite-density QCD with two flavors of Wilson 
quarks and the first calculation of the equation of state with 2 + 1 flavors of Wilson 
quarks. We are extending the studies towards our final objective of 2 + 1 flavor QCD 
directly at the physical point. 
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